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^  ABSTRACT 

The  ability  to  model  random  time  series  plays 
a  prominent  role  in  a  varlecy  of  applicationa  as 
exemplified  by  seismic  data  analysis,  dopplar 
radar  processing,  speech  processing,  adaptive 
filtering,  and,  array  processing.  Undoubtedly, 
cwo  of  che  more  popular  procedures  for  effecting 
such  time  series  models  are  che  classical  Fourier 
(MA)  approach  and  the  maximum  entropy  (AR)  method. 
In  this  paper,  a  cheoredcal  comparison  of  these 
contemporary  procedures  with  a  more  general  ARMA 
method  will  be  made.  It  will  be  demonstrated  Chat 
the  spectral  estimation  performance  of  the  ARMA 
method  typically  exceeds  that  of  its  more  special¬ 
ized  MA  and  AR  counterparts.  With  this  supremacy 
thus  established,  a  recently  developed  method  for 
effectively  generating  ARMA  model  estlmatea  from 
time  series  observations  will  be  then  presented. 

\ 

I .  INTRODUCTION 

A  problem  which  arises  in  a  variety  of  appli¬ 
cations  is  that  of  estimating  che  statistical 
characteristics  of  a  random  wlde-sense  stationary 
time  series  tx(n)i.  This  estimation  is  typically 
based  upon  a  finite  set  of  time  series  observations. 
In  many  signal  processing  applications,  only  the 
second  order  statistics  as  represented  by  the  time 
series'  autocorrelation  sequence 

rx(n)  *  Elx(n+m)x*(m) }  (1) 

is  required.  In  this  expression,  the  symbols  E 
and  *  denote  the  operations  of  expectation  and  com¬ 
plex  conjugation,  respectively.  Upon  taking  the 
Fourier  transform  of  this  deterministic  auto¬ 
correlation  sequence,  we  obtain  the  associated 
power  spectral  density  function 

S  (e^“)  *  l  r  (n)e  (2) 


It  often  happens  that  the  essential  attributes  of 
a  time  series  are  more  discernible  from  its  fre¬ 
quency  domain  spectral  density  function  than  from 
its  equivalent  time  domain  autocorrelation  sequence. 
It  is  with  this  in  mind  that  interest  in  spectral 
estimation  techniques  has  evolved. 
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In  the  classical  estimation  problem,  it  is 
desired  to  achieve  an  estimation  of  the  spectral 
density  function  (2)  from  observations  of  s  time 
series.  Without  loss  of  generality,  these  obser- 
vetions  will  be  here  taken  to  be  the  N  contiguous 
elements 

x(l),  x(2) . x(N)  (3) 

A  variety  of  procedures  have  been  proposed  for 
using  these  observations  to  achieve  a  spectral 
density  estimate.  Without  doubt,  the  overwhelming 
number  of  procedures  ultimately  result  in  a 
rational  spectral  density  model  which  fits  the  form 
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The  and  b^  coefficients  of  this  model  are 

referred  to  as  its  autoregressive  and  moving 
average  coefficients,  respectively.  This  model  is 
commonly  referred  to  as  an  autoregressive -moving 
average  (ARMA)  spectral  model  of  order  (p,q)  where 
q  and  p  denote  the  orders  of  the  numerator  and 
denominator  polynomials,  respectively.  It  is 
readily  shown  that  any  continuous  (in  u)  spectral 
density  may  be  approximated  arbitrarily  closely  by 
the  above  rational  model  If  the  order  (p,q)  is 
selected  suitably  large.  Thus,  che  robustness  of 
this  rational  model  is  apparent. 

In  studies  directed  cowards  spectral  analysis, 
the  preponderance  of  effort  has  been  directed  to¬ 
wards  tvo  special  cases  of  che  general  ARMA  model 
(A).  They  are  the  moving  average  (MA)  model  for  _ 
which  Ap(eJ“)  3  1,  and,  che  autoregressive  (AR)  p 
model  for  which  Bq(eJu).  *  bg.  The  spectral  density  — 
arising  from  a  MA  model  is  seen  to  contain  no  poles, 
and,  as  such  it  is  known  as  an  all-zero  model. 
Similarly,  the  AR  model  is  referred  to  as  an  all¬ 
pole  modal,  and,  the  general  ARMA  model  is  seen  to 
be  a  pole-zero  model.  Undoubtedly,  the  primary  TL 

reasons  for  Interest  in  che  special  case  MA  and  AR  _ 

models  are  chat  they: 


(1)  are  amenable  to  a  tractable  analysis, 
(il)  typically  provide  adequate  spectral 
•stlmatlon  performance. 
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(ill)  are  synthesizable  by  efficient  algorithm*. 

Despite  these  factors,  it  is  widely  recognized  that 
the  more  robust  ARMA  model  typically  provides  a 
much  superior  spectral  estimation  performance  and 
uses  fewer  model  parameters  in  the  estimate.  The 
main  impediment  to  its  use  on  a  wider  scale  has  been 
the  lack  of  a  specific  procedure  for  obtaining  the 
ARMA  model  parameters  in  a  computationally  effici¬ 
ent  manner.  Recently,  an  ARMA  spectral  modeling 
method  which  possesses  this  computational  capabili¬ 
ty  has  been  developed  [  1  ]  -  [  4 ] .  The  main 
features  of  this  new  procedure  are  outlined  in 
Sections  XV  and  VI. 

In  the  next  three  sections,  theoretical  pro¬ 
cedures  for  generating  MA,  AR,  and  ARMA  models  from 
a  finite  set  of  actual  autocorrelation  values  are 
presented.  This  is  then  followed  by  an  example 
section  which  treats  the  classical  problem  of 
spectral  estimating  a  time  series  composed  of  two 
sinusoids  in  additive  white  noise.  It  is  there 
shown  that  the  more  general  ARMA  modeling  procedure 
easily  outperforms  its  special  case  MA  and  AR 
modeling  procedures.  With  the  ARMA  spectral  raiallng 
methods  superior  performance  thereby  demonstrated 
for  this  idealistic  situation  (i.e.,  actual  auto¬ 
correlation  values  are  given) ,  we  next  direct  our 
attention  to  the  more  practical  problem  of  generat¬ 
ing  ARMA  spectral  estimates  from  a  finite  set  of 
observations  (3).  In  particular,  the  recently 
developed  high  performance  ARMA  modeling  method  is 
outlined  [2]&[4],  This  novel  procedure  has 
been  found  to  outperform  such  alternate  ARMA 
modeling  procedures  aa  the  Box-Jenklns  method [10] 
and  whitening  filter  approaches  [11 ] & [12] . 

II.  MA  SPECTRAL  MODELING 

There  exist  a  variety  of  procedures  for  obtain¬ 
ing  a  MA  model  of  a  wide  sense  stationary  time 
series.  These  Include  the  perlodogram  and  differ¬ 
ent  versions  of  the  autocorrelation  method.  In 
each  case,  the  spectral  estimate  will  be  of  the 
form 

Sx(eJu)  -  |bQ  +b1e_Ju,-f-  .  .  .+bqa-J<I“|2  (5) 

This  MA  spectral  model  requires  a  rather  large 
value  for  the  order  parameter  q  to  enable  it  to 
achieve  a  desirable  frequency  resolution  perform¬ 
ance.  Unfortunately,  this  requirement  can  lead  to 
a  rather  poor  spectral  estimation  behavior  in  the 
case  of  moderate  length  time  series  observations. 

This  behavior  will  be  illustrated  in  the  numerical 
example  section. 


in  which  w(n)  is  a  symmetric  window  function  which 
is  selected  to  effect  some  desired  behavior.  In 
the  pure  truncated  case,  the  rectangular  window  is 
used  in  which  case  w(n)  -  1  for  |q|<n.  Since  the 
autocorrelation  sequence  is  a  complex  conjugate 
symmetric  function  of  a  (i.e.,  rx(-n)  -  rx*(n)), 
one  can  readily  show  that  expression  (6)  can  be 
equivalently  represented  in  the  form  specified  by 
expression  (5)  whereby  the  prevailing  coefficients 
are  related  by 

q 

w(n)r  (n)  -  £  bvbk-n  for  °-nsq  (7> 

k-n 

Given  the  q+1  values  of  the  product  w(n)rx(n),  one 
may  readily  solve  this  nonlinear  system  of  q  +  1 
equations  to  obtain  the  MA  spectral  models  bfc  co¬ 
efficients  as  used  in  expression  (5).  Expressions 
(6)  and  (7)  then  constitute  a  systematic  procedure 
for  generating  a  MA  modal  of  time  series  baaed  on 
given  autocorrelation  values. 

In  most  applications,  one  has  available  only 
a  sat  of  time  series  observations  (3)  (and  not 
autocorrelation  values)  upon  which  to  generate  a 
MA  spectral  model.  If  expression  (6)  is  to  be  used 
for  this  objective,  it  will  then  be  necessary  to 
obtain  estimates  of  the  autocorrelation  elements 
from  the  given  time  series  observations.  The  un¬ 
biased  estimator  as  specified  by 

i  N~n 

rxC»)  *  jjlj  l  x(n+k)xTk)  0m<q  (8) 


la  often  used  for  this  purpose  In  which  It  Is 
assumed  that  the  MA  model  order  la  such  that  q<N. 
Alternatively,  the  perlodogram  method  has  served 
as  a  useful  procedure  for  effecting  a  MA  spectral 
estimate  from  a  set  of  contiguous  time  series  obser¬ 
vations  [  6]  .  The  perlodogram  possesses  the 
additional  advantage  of  being  efficiently  imple¬ 
mented  by  the  fast  Fourier  transform  algorithm. 

The  inherent  order  of  a  MA  perlodogram  model  is  N-l. 

III.  AR  SPECTRAL  MODELING 


In  this  section,  the  method  of  linear  pre¬ 
diction  Is  uaed  for  generating  an  AR  spectral  model 
associated  with  a  given  time  series.  In  particular, 
the  coefficients  of  the  pch  order  AR  spectral  model 
as  specified  Ipy 


y.*> 
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In  this  section,  we  will  be  concerned  with 
evolving  a  MA  spectral  estisate  procedure  for  the 
special  case  in  which  one  has  available  the  q+1 
autocorrelation  elements  rx(0) ,  rx(l),.  •  •  »rx(q). 

With  this  information  provided,  a  atandard  esti¬ 
mation  ia  generated  by  the  truncated  series  [5]. 

S  (eJ“)  -  l  w(n)r  (n)e'3“n  (6) 

*  n— q 


will  be  determined  by  solving  a  system  of  p+1 
linear  equations  in  the  p+1  coefficient  unknowns 
•1>  *2»  •  •  •,  ap,  bg.  These  equations  are  obtain¬ 
ed  by  considering  the  specific  problem  of  pre¬ 
dicting  the  time  series  element  x(n)  by  a  linear 
combination  of  the  p  most  recant  tins  sarlas 

alaaanta  x(n-l) ,x(n-2) . x(n-p).  It  will 

turn  out  chat  the  resultant  sat  of  optimal  equat¬ 
ions  thus  obtained  will  correspond  exactly  with 
those  aquations  which  arise  when  using  the 


2 


"maxinum  entropy”  method  of  spectral  estimation. 
This  equivalency  has  been  previously  recognized  [  7 ]. 

As  Indicated,  the  task  at  hand  Is  that  of 
estimating  the  time  series  element  x(n)  by  means  of 
che  linear  combination 

P 

x(n)  »  -  £  a^xCn-k)  (10) 

k*l 

In  which  che  optimal  prediction  coefficients  a^9 
are  to  be  ultimately  used  in  the  A R  spectral  model 
(9) .  The  prediction  error  Is  formally  given  by 


e(n)  -  x(n)  -  x(n) 

P 

-  x(n)  +  £  aJtx(n-k)  (11) 

k-1 

Ic  Is  now  desired  to  select  the  complex  valued  a^ 
prediction  coefficients  so  chat  che  expected  value 
of  the  error's  magnitude  squared  is  minimized.  One 
may  straightforwardly  show  thac  this  mean  squared 
error  is  given  by 

,  P 

E{  i  e(n)  |  ’  -rx(0)  +J  (Vx('k)  *Vx(k)l 
P  P 

+  l.  Uvx("-kl  (12) 

k-i  m-1 

Upon  using  standard  calculus  methods,  it  Is  found 
that  the  minimizing  predictor  coefficients  must 
satisfy  the  following  system  of  p  linear  equations 

P 

1  a."  r  (m-k)  »  -r  (m)  for  l<raip  (13) 
k-l  x 


If  this  optimal  selection  is  inserted  into  ex¬ 
pression  (12),  the  minimum  mean  squared  error  is 
found  to  be 

,  P 

Ei  ;  e9  (n)  |  j  •  r  (0)  +  £  r  (-k)  (14) 

x  k-1  x 


In  suosary,  the  optlisal  predictor  coefficients 
are  obtained  upon  solving  che  system  of  linear 
equations  (13),  and,  Che  corresponding  minimal 
mean  squared  error  is  computed  by  means  of  relation¬ 
ship  (14).  From  a  computational  viewpoint,  however, 
a  more  efficient  method  for  iteratively  obtaining 
the  optimal  predictor  performance  is  obtained  by 
incorporating  relationships  (13)  and  (14)  into  the 
single  expression 


;rx(0)  rx(-l) 
r  (1)  r  (0) 

X  X 


rx(-p) 

rx(-p+l) 


>1 


(15) 


Vp) 


rx(p-l)....  rx(0) 


0 


in  Which  Ep  -  E{ [ •* <n) | 2 > .  One  may  solve  this  Toeplitz 
system  of  equations  using  the  Levinson  algorithm 
in  a  computationally  efficient  order  update  manner 


(e.g.,  see  ref.  [7]).  Namely,  che  optimal  co¬ 
efficients  of  the  p+l8t  order  predictor  may  be 
recursively  obtained  from  the  optimal  coefficients 
of  the  pth  order  predictor.  As  indicated  previous¬ 
ly,  the  system  of  equations  (15)  is  identical  to 
that  obtained  when  using  che  maximum  entropy 
method  of  spectral  estimation. 

If  che  optimal  predictor  is  performing  its 
objective,  it  follows  that  the  prediction  element 
x(n)  will  contain  ail  which  is  predictable  in  x(n) . 
As  such,  the  prediction  error  (11)  is  white  like 
in  behavior  and  its  spectral  density  is  then  given 
by  Se(eJ“)  -  Ep.  This  behavior  is  of  course  de¬ 
pendent  on  making  the  order  parameter  p  sufficient¬ 
ly  large  so  as  to  achieve  the  desired  perfect 
prediction.  Assuming  this  prediction  behavior,  it 
then  follows  from  relationship  (11)  that  the 
spectral  density  function  of  the  time  series  (x(n)) 
is  given  by  expression  (9)  in  which  the  coefficient 
bo2  -  Ep  and  the  a^  coefficients  are  obtained 
upon  solving  relationship  (13)  or  equivalently 
relationship  (15). 

In  this  analysis,  it  has  been  assumed  that 
one  has  available  the  time  series  autocorrelation 
values  rx(0),  rx(l) ,  ....  rx(p).  More  realistic¬ 
ally,  such  perfect  autocorrelation  knowledge  is 
almost  never  available.  In  a  typical  application, 
one  has  available  only  a  sampled  set  of  time  series 
observations  as  exemplified  by  expression  (3).  If 
the  AR  spectral  estimation  procedure  as  represented 
by  expression  (15)  is  to  be  incorporated,  one  could 
use  the  given  time  series  observations  to  obtain 
estimates  of  che  required  autocorrelation  elements. 
Alternatively,  a  set  of  deterministic  prediction 
error  equations  can  be  minimized  so  as  to  obtain 
the  prediction  coefficients.  In  effect,  this  is 
the  approach  usually  taken  in  evolving  the  Burg 
algorithm  [ 7  ) . 


IV.  ARMA  SPECTRAL  MODELING 

The  time  series  (x(n))  is  said  to  be  an  ARMA 
process  of  order  (p,q)  if  ic  is  generated  according 
to  the  linear  causal  relationship 

p  q 

x(n)  +  £  a.  x(n-k)  -  £  b.w(n-k)  (16) 

k-1  k-0 

where  the  excitation  sequence  (w(n) }  is  a  zero  mean 
white  noise  time  series  whose  individual  elements 
have  variance  one.  It  is  a  simple  matter  to  show 
thac  che  spectral  density  function  associated  with 
this  response  time  series  (x(n)}  is  given  precisely 
by  expression  (4).  Thus,  there  is  an  equivalency 
between  a  rational  spectral  density  model  and  the 
response  of  a  linear  system  to  a  white  noise 
excitation. 

a^  Coefficient  Determination 

A  procedure  for  Identifying  an  ARMA  model's 
ak  autoregressive  coefficients  involves  examining 
its  second  order  statistical  characterization. 

This  is  achieved  by  first  multiplying  both  sides  of 
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expression  (16)  by  x*(n-m)  and  then  taking  expect¬ 
ed  values.  The  results  of  this  operation  are  the 
well  known  Yule-Walker  equations 

P 

r  (m)  +  l  a.  r  (m-k)  *  0  for  ts  >  q  (17) 
x  fc-1 

where  It  Is  important  to  note  that  the  lag  variable 
a  Is  here  restricted  to  be  larger  than  the  ARMA 
model's  numerator  order  q. 1 

In  effect,  the  above  Yule-Walker  equations 
Indicate  that  an  ARMA  time  series'  autocorrelation 
elements  are  Interrelated  In  a  well  defined  linear 
manner  for  appropriate  lag  values.  This  obser¬ 
vation  then  provides  the  vehicle  for  determining 
che  A MW  Jdel's  associated  a^  coefficients.  To 
be  more  sveclflc,  let  us  now  express  che  first  t 
Yule-Walker  equations  (17)  In  che  following  matrix 
format 


rx(q) 

rx(q-D.  • 

'  -rx(9+l-p) 

al 

rx(q+D 

,rx(q+2) 

rx(q)  •  • 

.  ,rx(q+-2-p) 

a2 

m  — 

rx(q+2) 

SP 

rx(q+t-l) 

rx(q+t-2) . 

•  -rx(9+t-P) 

v  (q+e) 

* 

(18) 

or  In  the  more  coapact  representation 

R  a  -  -r  (19) 

where  R  Is  a  t»p  matrix  and  a^  and  £  are  pxl  and 
t*l  vectors,  respectively. 

It  Is  readily  shown  chat  this  system  of  aquations 
will  have  a  unique  solution  provided  that  t  £  p 
To  obtain  che  ARMA  model's  afc  coefficients,  one 
then  slaply  solves  this  system  of  linear  equations. 
From  a  computational  viewpoint.  It  Is  convenient 
to  set  c  equal  to  Its  minimum  value  of  p.  For 
reasons  which  will  be  spelled  out  In  Section  VI, 
however.  It  will  often  be  advantageous  to  let 
t  cake  on  values  exceeding  p. 

b^  Coefficient  Determination 

To  determine  the  ARMA  model's  b^  coeffici¬ 
ents,  It  will  be  expedient  to  Introduce  che  auto¬ 
correlation  sequence's  causal  Image 

r^"(n)  •  rx(n)u(n)  -  ^  rx(0)«(n)  (20) 

where  u(n)  and  S(n)  denote  the  etandard  unlt- 
etep  and  Kroneckar  delta  sequences,  respectively. 
The  autocorrelation  sequence  may  be  recovered  from 
its  causal  Image  through  the  relationship 

rx(n)  -  rx+(n)  +  rx+(-n)*  (21) 

*The  Yule-Walker  equations  associated  with  lag 
valuea  0£m£q  will  Involve  the  ARMA  model's  bfc 
coefficients'  In  a  nonlinear  manner. 


whose  validity  Is  established  by  making  use  of  che 
complex  conjugate  syiwetry  property  of  auto¬ 
correlation  sequences  (i.e.,  r(-n)  -  r*(n)).  The 
Fourier  transform  of  expression  (21)  Is  found  to 

be 

Sx(eJ“)  -  S*(e3u)  +  Sx+(eJ")  (22) 

in  which  Sj^(e^“)  denotes  the  Fourier  transform  of 
the  autocorrelations'  causal  Image  (20).  In  what 
is  to  follow,  a  systematic  procedure  for  Identify¬ 
ing  S±(eJ“)  will  be  given  which, with  the  utilisat¬ 
ion  of  expression  (22),  results  In  the  overall  times 
series'  spectral  density. 


This  Identification  Is  perhaps  best  achieved 
by  Introducing  the  following  auxiliary  sequence 

P 

c(n)  *  r +(n)  +  T  a.r  (n-k)  (23) 

*  k-1  x 

in  which  the  a^  coefficients  as  generated  from 
expression  (19)  are  used.  Due  to  the  nature  of  the 
causal  image  sequence  and  the  underlying  Yule- 
Walker  equations  (17) ,  it  Is  readily  shown  that 
this  auxiliary  sequence  is  identically  zero  outside 
the  indexing  range  0  £  n  £  max(p.q).  Using  this 
fact,  the  Fourier  transform  of  relationship  (23) 
is  found  to  yield 
s 

C  (e^“)  -  l  c(n)e”^un  ,  s  -  max(p.q)  (24) 
®  n-0 

-  A(.JXVU)  (25) 

P  * 

Using  this  result  in  aquation  (22),  the  required 
ARM*  spectral  density  formulation  la  obtained 


.  A  (el“)C*(.Ju)  -aV“)C  (eJ“) 
<5  ( ,_2 - ! - _JE - ! - 

x  A(e^)A*(.J“) 

P  P 


To  obtain  the  ARMA  model's  bfc  coefficients, 
we  next  incorporate  relationships  (4)  and  (26)  to 
generate  the  complex  conjugate  sysaaetrical  poly¬ 
nomial  (in  the  eJ“)  expression 


B  (eJ")B*(«j“)  -  A  (e>)C*(«J“)  ♦  A*(.J")C,(e>) 

a  q  ps  F  * 


(27) 


A  spectral  factorazstlon  of  the  right  side  poly¬ 
nomial  will  yield  2q  roots  which  occur  In  complex 
conjugate  reciprocal  sets.  One  then  need  only 
select  an  appropriate  q  of  these  roots  to  deter¬ 
mine  the  required  Bq(eJ“)  term  (e.g.,  the  minimum 
phase  selection). 


To  summarize,  the  spectral  density  function 
and  the  model  coefficients  corresponding  to  an 
ARMA  time  series  of  order  (p,q)  may  be  obtained  by 
following  the  systematic  procedure  outlined  in 
Table  1.  To  carry  out  this  process,  it  is 
necessary  to  have  knowledge  of  the  order  pair 
(p,q)  and  the  autocorrelation  elements  rx(n)  for 
0  <  n  <  q+p. 
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1. 

Solve  Relationship  (19)  for  the  p  auto¬ 
regressive  ag  coefficients.  This  will  re¬ 
quire  setting  t  _>  p. 

2. 

Cenarate  the  auxiliary  sequence  c(n)  and  its 
Fourier  transform  using  expressions  (23)  and 

(24),  respectively. 

3. 

The  required  spectral  density  is  glvan  by 
relationship  (26). 

4. 

Perform  a  spectral  factorization  of  the 
polynomial  B(eJo»)B*(eJ“)  as  given  by  ex¬ 
pression  (27)  to  obtain  the  required  bg 
coefficients. 

Table  1:  Generation  of  the  spectral  density  and 
che  ARMA  model  parameters  associated  with  a 
given  set  of  autocorrelation  valuea. 

V.  EXAMPLE 

To  demonstrate  the  relative  effectiveness  of 
che  MA,  AR,  and  ARMA  modeling  schemes  presented  in 
che  last  chree  sections,  let  us  consider  the  class¬ 
ical  problem  of  generating  spectral  estimates  of 
a  time  series  composed  of  two  sinusoids  in  additive 
whlce  noise.  Namely,  che  time  series  will  be 
governed  by  che  relationship 

x(n)  »  a^sin(u^n+0^)  +  a^slnlu^n+ej)  +  w(n)  (28) 

in  which  6}  and  are  Independent  random  varlablea 
uniformly  distributed  on  [0,211]  and  v(n)  is  a  zero 
mean  white  noise  process  whose  individual  elements 
have  variance  a  .  The  sinusoidal  amplitudes  a^  and 
a2,  and,  normalized  frequencies  and  uj  will  be 
here  taken  to  be  unknown  constants.  2 


It  is  readily  shown  chat  Che  autocorrelation 
sequence  corresponding  to  this  time  series  is 
given  by 

,2  .2  2 

r  (n)  »  1  cos(u.a)  ♦  2  coa(u,n)  +  0*6(0)  (29) 

*  T  1  T  2 

Upon  caking  the  Fourier  transform  of  this  auto¬ 
correlation  sequence,  che  associated  spectral  den¬ 
sity  tunc  cion  is  found  to  be 
.  ,  a  2 

Sx(eJ<“)  »_l_n[4(»i-(i>^)+4(ivHi^)  ] 

2  a  2  2 

+  2  fl  f  6  (is— )+6  (uy+to..)  T  +0 

2  2  i 

for  -II  <.n  (30) 

This  spectral  density  function  is  seen  to  be  com¬ 
posed  of  dlrac  delta  functions  located  at  fre¬ 
quencies  t  m  4  i  wj  riding  on  top  of  a  constant 
value  t2  due  Co  che  additive  white  noise. 


choices  for  the  time  series  parameters  were  taken 
to  be  aj  “  /2Q,  wj  m  0.4H,  S2  “  ^2,  U2  *  0.426II  and 
o2  m  1.  This  selection  provides  individual  sinu¬ 
soid  signal- co-noise  ratios  of  lOdB  (decibels)  and 
OdB.  Due  to  the  relative  closeness  of  the  sinu¬ 
soid  frequencies,  this  example  provides  an  excel¬ 
lent  measure  of  the  frequency  resolution  capabili¬ 
ties  of  the  three  modeling  procedures.  A  brief 
description  of  the  results  obtained  for  this 
example  now  follows. 

MA  Spectral  Estimates 

The  autocorrelation  elements  as  specified  by 
expreaslon  (29)  where  Incorporated  into  the  MA 
spectral  model  relationship  (6)  in  which  the 
window  function  is  caken  to  be  rectangular  (l.e. , 
w(n)  ■  1  for  0  £  n  £  q) .  The  spectral  estimates 
thereby  achieved  for  che  specific  order  selections 
q  «  15,  30,  and  200  are  displayed  in  Figure  1. 

From  these  plota  it  is  clear  that  the  MA  modeling 
procedure  ia  unable  to  resolve  che  sinusoids  for 
che  order  selections  q  *  15  and  30.  Whan  q  la  set 
to  200,  it  is  possible  to  Just  barely  detect  the 
presence  of  the  lower  amplitude  sinusoid.  It  Is 
apparent  from  this  example  that  classical  Fourier 
approaches  provide  relatively  poor  vehicles  for 
achieving  frequency  resolution  even  when  exact 
autocorrelation  elements  are  used. 

AR  Spectral  Estimates 

The  autocorrelation  elements  (29)  were  next 
Incorporated  into  the  optimum  one-step  predictor 
(or  msxlmun  entropy)  expression  (IS) ,  with  AR 
order  choices  of  p  -  15  and  30.  The  two  AR  spect¬ 
ral  estimates  which  resulted  are  shown  in  Figure 
2  where  It  is  apparent  chat  a  frequency  resolution 
is  achieved  for  p  -  30,  but,  not  for  p  -  IS.  In 
contrast  to  the  classical  Fourier  approach,  the 
AR  spectral  modeling  procedure  is  capable  of 
achieving  the  required  frequency  resolution  with 
a  reasonably  small  order  model. 

ARMA  Spectral  Estimate 

In  the  final  modeling  approach,  the  auto¬ 
correlation  elements  (29)  were  next  used  to  obtain 
an  ARMA  model  of  order  (4,4)  using  the  procedure 
as  outlined  in  Table  1.  The  resultant  spectral 
estimate  is  plotted  in  Figure  3  and  is  seen  to 
correspond  precisely  with  that  as  given  in 
equation  (30). 3  This  should  not  be  surprising 
considering  che  fact  that  the  given  autocorrelation 
sequence  (29)  is  an  ARMA  time  aeries  of  order  (4,4). 
Thus,  the  model  used  in  this  case  precisely  matches 
the  time  series  being  examined. 

VI.  HIGH  PERFORMANCE  METHOD  OF 
ARMA  SPECTRAL  MODELING 


Using  che  given  autocorrelation  elements  (29) 
as  entries,  the  three  spectral  estimation  pro¬ 
cedures  just  described  were  next  utilized  to 
generate  MA.  AR,  and  ARMA  models.  The  specific 


2The  normalized  frequency  variables  era  taken  to 
lie  In  Che  interval  (0 , ft ] . 


It  is  possible  to  adapt  many  of  the  ideas  of 
Section  IV  to  achieve  an  ARMA  spectral  estimate 
whan  only  chs  time  series  observation  (3)  are 


If  infinite  precision  computatlonswere  utilised, 
thara  would  be  an  exact  pole-zero  cancellation  and 
an  associated  constant  spectral  density  function. 
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available  (and  not  autocorrelation  value*).  We 
shall  again  treat  separately  the  cases  of  auto¬ 
correlation  and  moving  average  coefficient  deter¬ 
mination. 

Autoregressive  Coefficient  Estimation 


ARMA  model  order  choice.  In  any  case,  the  system 
of  equations  with  these  estimate  substitutions 
will  give  rise  to  the  t«l  Yule-Walker  approximation 
error  vector  as  specified  by 

e  •  Y+X  a  +  Y+x  (36) 


To  implement  the  autoregressive  coefficient 
selection  process  as  represented  by  relationship 
(19) ,  1c  will  be  necessary  to  compute  appropriate 
autocorrelation  estimates  from  the  given  sat  of 
clme  series'  observations.  The  high  performance 
ARMA  method  effects  these  estimates  in  Che  guise 
of  a  convenient  matrix  format  which  lends  itself 
to  a  particularly  efficient  computational  realiza¬ 
tion  [2)&[4],  In  particular,  Che  autocorrelation 
matrix  and  vector  required  in  expression  (19)  are 
estimated  according  Co 

R  *  YfX 

(31) 

1  * 

(32) 

where  the  dagger  symbol  +  denotes  the  operation  of 
complex  conjugate  transposition.  The  (N-p)«p 
Toeplitz  type  matrix  X  is  specified  by 

'x(pl  x(p-l).,  .  . 

x(l) 

x(p+l)  X (p )  .  .  . 

X  -  1  . 

x(2)  I 

1 

.  i 

1 

(33) 

1  x  ( N-  L )  x(N-2) 

.  I 

x(N-p)  | 

J 

while  the  (N-p)*c  Toeplitz 
form 

r"x(p-q)  x(p-q-l)  . 

type  matrix  Y 

.  .x(p-q-t+l)"j 

has  the 

x(p-q+l)  x(p-q)  .  . 

.  • x(p-q-t+2) 1 

I 

Y  * ! 

(34) 

jx(N-q-l)  x(N-q-2)  . 

.  . x(N-q-t) 

and  x  is  a  (N-p)«l  vector 

given  by4 

x  *  |x(p+l) ,  x(p4-2) , 

.  .  .  x(N) ] ’ 

(35) 

In  formulating  matrix  Y,  we  have  used  the  con¬ 
vention  of  setting  to  zero  any  elements  x(k)  for 
which  k  lies  outside  the  observation  index  range 
1  <  k  «  N. 

If  the  autocorrelation  matrix  and  vector 
estimates  (31)  and  (32),  respectively,  are  substi¬ 
tuted  into  the  Yule-Walker  relationship  (19), 
however,  it  is  generally  found  that  the  resultant 
system  of  t  equations  in  the  p  autoregressive 
coefficients  Is  inconsistent  for  t  >  p.  This  is 
due  to  inevitable  inaccuracies  in  the  auto¬ 
correlation  estimates,  and,  to  a  possible  Improper 
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A  more  generalized  version  of  this  estimation 
schema  can  be  obtained  by  substituting  the  Inte¬ 
ger  k  for  p  wherever  p  appears  in  relationships 
(33)  -  (35).  For  ease  of  presentation,  k  is 

here  restricted  to  be  p. 


Upon  caking  the  expected  value  of  5,  It  Is  found 
that  for  the  ARMA  modeling  order  choice  In  which 
q  >.  p,  that  this  expectation  results  in 


E(e(k) } 


(N-q 


1-  P 

-k)  rx(q+k)  +  [  a^fq+k- 


m“l 


1  <_k£t 
(37) 


while  for  the  modeling  order  case  q  <  p  this  ex¬ 
pectation  produces 


E(e(k)}  - 


r  p 

r  (q+k)  +  J  a  r 


(N-p)  jr  (q+k)  +  [  a  r  (<^k-m) 


m-l 


(N-q-k) ir  (q+k)  +  la  r  (q+k-m) 


,  l£k£p-q 

,  p-q  <  k  <_t 

(38) 


In  either  ordering  case.  It  is  seen  that  when  the 
time  series  is  an  ARMA  process  of  order  (p,q),  the 
expected  value  of  the  error  vector  e  can  be  made 
equal  to  zero  by  a  proper  choice  of  the  auto¬ 
regressive  coefficient  vector  a  .  Namely,  this 
selection  would  be  such  that  the  underlying  Yule- 
Walker  equations  (17)  are  satisfied. 5  This  implies 
that  the  system  of  equations  (36)  with  e  -  £  pro¬ 
vides  an  unbiased  and  a  consistent  estimate  of  the 
Yule-Walker  equations. 


With  the  above  thoughts  In  mind,  an  appealing 
approach  to  selecting  the  autoregressive  co¬ 
efficient  vector  Is  lsMdlately  suggested.  Namely, 
a  is  chosen  so  as  to  make  the  error  vector  "as 
close”  to  Its  expected  value  of  £  as  possible. 

This  is  of  course  predicated  on  the  assumption 
that  the  time  series  is  an  ARMA  process  of  order 
(p,q)  or  less.  In  order  to  attain  a  tractable 
procedure  for  selecting  an  appropriate  auto¬ 
regressive  coefficient  vector,  we  shall  Introduce 
the  following  quadratic  functional 

f (a)  -  *+A  *  (39) 


in  which  A  is  a  t*t  positive-definite  diagonal 
matrix  with  diagonal  elements  1^  Chat  is 
Introduced  In  order  to  provide  one  with  Che  option 
of  weighting  differently  the  various  error  vector 
components.  It  is  a  simple  matter  to  show  that  an 
autoregressive  coefficient  vector  which  will 
render  this  quadratic  functional  a  minimum  must 
satisfy 

XfY  AY+X  a*  -  -  X+Y  AY+x  (40) 


J -  , 

A  little  though  will  convince  oneself  that  this 

same  conclusion  will  be  reached  if  both  q  and  p  j 

arc  at  least  equal  co  the  numerator  and  danoml-  j 

nator  orders,  respectively,  of  the  underlying  j 

ARMA  time  serlas.  I 
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On*  than  sis ply  solve*  this  consistent  system  of 
p  linear  equations  In  the  p  unknown  autoregressive 
coefficients  to  obtain  an  estlaat*  for  the  denomi¬ 
nator  dynamics  of  the  IRMA  model. 

Moving  -  Average  Coefficient  Estimation 

Thera  exist  several  procedures  for  estimating 
the  ARMA  nodal '*  moving  average  coefficients.  Ua 
shall  now  briefly  describe  two  procedures  which 
have  provided  satisfactory  performance  and  in  a 
sense  complement  one  another. 

(1)  Method 

The  procedure  which  has  provided  che  beat 
frequency  resolution  behavior  is  a  direct 
adaption  of  the  method  as  described  in 
Section  IV  (and  ref.  [3])-  In  particular, 
using  che  set  of  autoregressive  coefficient 
estimates  as  obtained  from  expression  (40) 
and  a  suitable  set  of  autocorrelation 
estimates  rx(n)  for  n » 0,1, . . . ,mex(q,p)  one 
computes  the  c(n)  coefficients  using  expres¬ 
sion  (23) .  These  coefficients  are  than  used 
to  achieve  the  desired  ARMA  spectral  estimate 
whan  Incorporated  Into  relationship  (24)  and 
ultimately  relationship  (26).  Although  pro¬ 
viding  an  excellent  frequency  resolution 
behavior,  this  procedure  suffers  the  drawback 
of  not  having  a  guaranteed  nonnegative  defin¬ 
ite  spectral  density  function.6  It  la  with 
this  In  mind  that  the  following  wall  known 
smoothed  perlodogram  method  was  adapted  [13]. 


ek(n)  *  £(n+p  +  l+kd)  0£n<^q  (43) 

0  <k  <L-1 

where  "d"  is  a  positive  Integer  which  specifies  the 
time  shift  between  adjacent  segments.  These  Indi¬ 
vidual  segments  will  overlap  if  d  and  will  per¬ 
fectly  partition  the  residual  sequence  when  d  “  q  + 1. 
In  order  to  Include  only  computed  elements,  the 
relevant  parameters  must  be  selected  so  that 
q  +  p  +1  +  (L-l)d  <^H.  Next  the  perlodogram  for  each 
of  these  L  segments  Is  taken  and  then  averaged  to 
obtain  the  desired  qch  order  smoothed  perlodogram, 
that  is 

s  (eJ“)  Y|=rl  r  w(n)cu(n)e-J““l  \  (44) 

E  L  k-0 


q 

2 

1 

q+1 

l  w(n)t.  (n)e_^"° 
n«0 

where  w(n)  la  a  window  sequence  that  Is  normally 
selected  to  be  rectangular  (l.e.,  w(n)  * l//q+l  for 
0  <n  £q) .  It  is  readily  shown  that  this  procedure 
results  In  a  desired  nonnegative  qcl>  order  MA 
spectral  density  estimate.  Unfortunately,  Its 
frequency  resolution  capability  is  generally  not 
of  the  same  quality  as  that  of  the  cfc  method.7 
On  the  other  hand,  the  smoothed  perlodogram  method 
provides  more  smoothly  behaved  spectral  estimates 
which  contain  fewer  spurious  effects. 

To  summarize,  the  required  ARMA  spectral  model 
is  obtained  by  following  the  systematic  procedure 
outlined  In  Table  2.  The  numerator  dynamic  esti¬ 
mation  procedure  to  be  used  will  of  course  depend 
on  the  particular  characteristic  being  sought  (e.g., 
frequency  resolution,  smoothness,  etc.) 


(ii)  Smoothed  Perlodogram  Method 


In  the  smoothed  perlodogram  approach, 
one  first  computes  che  so-called  "residual" 
time-series  elements  according  to  the  relation¬ 
ship  (see  ref.  (3]). 

P  „  . 

t(n)  -x(n)  +  J  a.  x(n-k)  for  p<  n<  N  (41) 
k-1 


In  which  the  a^*  autoregressive  coefficients 
as  obtained  by  solving  expression  (40)  are 
Incorporated.  From  this  relationship  the 
spectral  density  expression  directly  follows 


sK(.J“) 


St(.J“) 


(42) 


If  S„(*J“)  Is  to  correspond  to  an  ARMA  spect¬ 
ral  model  of  order  (p,q),  it  Is  clear  that  a 
qth  order  MA  spectral  estimate  for  the  resi¬ 
dual  spectral  density  St(eJ“)  must  be  obtained 
and  then  substituted  into  relationship  (42). 
The  smoothed  perlodogram  has  been  found  to  be 
a  useful  tool  for  this  purpose. 


lo  the  smoothed  perlodogram  method,  one  first 
partitions  the  computed  residual  el amenta  (41)  into 
L  segment*  each  of  length  q+1  as  specified  by 


i. 

Specify  values  for  the  ARMA  model's  order 
parameter  pair  (q,p),  the  Yule-Walker 
equation  parameter  t,  and,  the  weighting 
metrics' s  diagonal  elements 

2. 

Using  the  time  series  observations  x(l), 
x(2) , . . . ,x(N) ,  construct  the  matrices  X,  Y, 
and  vector  x  according  to  relationships  (33), 
(34),  and  (35),  respectively. 

3. 

Determining  the  model's  autoregressive  co¬ 
efficients  by  solving  relationship  (40). 

4. 

The  numerator's  dynamics  are  obtained  by 
using  either  the  (1)  cj.  method,  or, 

(11)  the  smoothed  perlodogram  method. 

Table  2:  Basic  steps  of  the  standard  high  per¬ 
formance  ARMA  spectral  estimation  method : 

The  Block  Processing  Mod*. 

The  improved  spectral  estimation  performance 
obtained  In  using  this  high  performance  method 
over  contemporary  ARMA  techniques  such  as  the  Box- 
Jenklna  method  Is,  to  a  large  extent,  a  conse¬ 
quence  of  selecting  the  integer  t  to  be  larger 
than  the  minimal  value  p.  With  the  corresponding 
larger  set  of  Yule-Walker  equations  that  are 


This  shortcoming  may  be  superficially  avoided  by 
taking  the  absolute  value  of  the  spectral  estlaat*. 


7 A  similar  approach  shares  tbs  same  attributes  as 
does  the  smoothed  perlodogram  [14], 
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thereby  being  approximated,  it  intuitively  follows 
that  the  model's  autoregressive  coefficients  will 
be  less  sensitive  to  autocorrelation  estimate 
errors  which  are  embodied  in  Y+X  and  Y+x  chan  would 
be  the  case  if  t  were  sec  to  p  (as  in  Che  Box- 
Jenklns  method) .  This  anticipated  improvement  in 
spectral  estimation  behavior  when  using  Che  high 
performance  method  has  in  fact  been  realized  on  a 
rather  large  number  of  numerical  examples  [1]-{4J. 

It  is  shown  in  reference  [4]  &  [8j  that  this  high 
performance  method  also  lends  Itself  to  a  particu¬ 
lar  fast  adaptive  implementation  mode  when  t  *  p . 
With  the  two  attributes  of  improved  spectral 
estimation  performance  and  computational  efficiency, 
this  new  procedure  promises  to  be  an  important 
spectral  estimation  tool. 

It  is  of  interest  to  note  that  when  q  -  0  and 
t  «  p,  the  high  performance  ARMA  spectral  esti¬ 
mation  method  reduces  to  the  well  known  AR  co- 
variance  method.  Moreover,  upon  letting  t  exceed 
p,  Che  resultant  set  of  expanded  AX  Yule-Walker 
equation  approximations  will  cypically  result  in 
better  spectral  estimates  than  the  standard  AX  co- 
variance  method,  to  the  author^  knowledge,  this  ap¬ 
proach  has  not  been  used  in  che  various  AX  spectral 
estimation  procedures  developed  to  date. 

VII.  CONCLUSION 

A  computationally  efficient  closed  form 
method  of  ARMA  spectral  estimation  has  been  pre¬ 
sented.  It  is  predicated  on  the  approximation  of 
a  set  of  Yule-Walker  equation  estimates  which  are 
generated  from  a  given  set  of  time  series  obser¬ 
vations-  The  ARMA  model's  autoregressive 
coefficients  are  determined  by  solving  a  consistent 
system  of  linear  equations. 

The  spectral  estimation  performance  of  this 
ARMA  modeling  procedure  has  been  empirically  found 
to  exceed  that  of  such  counterparts  as  Che  maximum 
entropy  and  Box-Jenklna  methods  (a.g.,  see  refs. 
[l]-[4]  &  [9]).  This  behavior  is,  to  a  large 
extant,  a  consequence  of  the  fact  that  more  than 
the  minimal  number  of  Yule-Walker  equation  esti¬ 
mates  are  being  approximated  to  obtain  the  result¬ 
ant  ARMA  model  parameters. 
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